## Replication code for Figure 1
##
## Christopher Adolph, Kenya Amano, Bree Bang-Jensen, 
## Nancy Fullman, Beatrice Magistro, Grace Reinke, Rachel Castellano,## Megan Erickson, and John Wilkerson, "The Pandemic Policy U-Turn:## The role of partisanship, public health, and race in decisions
## to ease COVID-19 social distancing policies in the U.S.", 
## forthcoming in Perspectives on Politics.
##
##
## Christopher Adolph
## 18 June 2021     faculty.washington.edu/cadolph

## Clear memory
rm(list=ls())


######################################################
####  Set up #########################################
## File name to save 
file <- "Figure1easingPrepolish"

## Use fancy fonts?  Requires setting the correct path below to GillSans font (on Macs)
fancy <- FALSE

# Set base directory
basedir <- getwd()
rowdata.folder <- basedir

######################################################

## Load libraries
library(tile)           # Graphics
library(RColorBrewer)   # Colors
library(lubridate)
library(readxl)
library(tidyverse)

## Set up colors
col <- brewer.pal(7, "Set1")
nicered <- col[1]
niceblue <- col[2]
nicegreen <- col[3]
nicepurple <- col[4]
niceorange <- col[5]
nicebrown <- col[7]


## Set some parameters affecting look of plots
xOffset <- 0.415 # original 0.415; .65
yOffset <- -3
cex <- 0.55
circleSize <- 1.15  # original 1.3
textSize  <- 0.59 # original 0.625

DateTypeList <- 
  c("DateEnacted",       
    "DateSubNonReligIndoorEasing")

StopDateType <-
  "DateSubNonReligIndoorEasing"


StatePolicyTable <-  
  read.csv(paste(rowdata.folder, "StatePolicy.csv", sep="/")) %>%   
  select(-SubstatePolicyDummy) %>% 
  filter(DateType %in% DateTypeList) %>% 
  arrange(StatePolicy, DateType, StateName)

#####################################################
##  Create wide table for plotting
##  Enacted date

StatePolicyEnactedTableWide <- 
  StatePolicyTable %>% 
  select(-SubstatePolicyDummy) %>% 
  filter(DateType == "DateEnacted") %>% 
  pivot_wider(names_from = StatePolicy, values_from = PolicyDate)

## Only count policy enact up to July 11

for (j in 4:8) {
  for (i in 1:nrow(StatePolicyEnactedTableWide)) {
    if (!is.na(StatePolicyEnactedTableWide[i,j]) && (StatePolicyEnactedTableWide[i,j]>20200711))
      StatePolicyEnactedTableWide[i,j] <- NA
  }
}

##  Policy eased date
StatePolicyStopTableWide <- 
  StatePolicyTable %>% 
  select(-SubstatePolicyDummy) %>% 
  filter(DateType == StopDateType) %>% 
  pivot_wider(names_from = StatePolicy, values_from = PolicyDate)


## Only count policy stops up to July 11

for (j in 4:8) {
  for (i in 1:nrow(StatePolicyStopTableWide)) {
    if (!is.na(StatePolicyStopTableWide[i,j]) && (StatePolicyStopTableWide[i,j]>20200711))
      StatePolicyStopTableWide[i,j] <- NA
  }
}


###################################################################
#Create function to identify a state that implement a policy but still haven't stopped
FindNotStopStateFunc <- 
  function(policy){
    StatePolicyTable %>% 
      filter(DateType %in% c("DateEnacted", StopDateType)) %>% 
      pivot_wider(names_from = DateType, values_from = PolicyDate) %>% 
      filter(is.na(DateEnacted)==F&is.na(!!sym(StopDateType))==T) %>% 
      filter(StatePolicy %in% policy) %>% 
      dplyr::select(StatePostal) %>% 
      as_vector()
  }



#Create function to identify a state that never implement a policy
FindNeverEnactStateFunc <- 
  function(policy){
    StatePolicyTable %>% 
      filter(StatePolicy %in% policy,
             DateType == "DateEnacted",
             is.na(PolicyDate) == T) %>% 
      dplyr::select(StatePostal) %>% 
      as_vector()
  }


## Weeks past function

weeksPast <- function(date, offset=4) {
  floor((date-offset)/7)
}



###########################################################################################
###########################################################################################

# PLOT

###########################################################################################

## First compute cumulative cases at risk by start of easing phase (April 6 2020)

## Attach data for convenience
attach(StatePolicyEnactedTableWide)

## Convert to days since 26 February 2020 (first reported US community transmission of COVID-19) 

# here I renamed stayathome and restaurantrestrict (they were previously shelterinplace and restaurantclose)
BarRestrictDays0 <- as.numeric(difftime(ymd(BarRestrict), ymd("20200226"), units="days"))
RestaurantRestrictDays0 <- as.numeric(difftime(ymd(RestaurantRestrict), ymd("20200226"), units="days"))
AnyBusinessCloseDays0 <- as.numeric(difftime(ymd(AnyBusinessClose), ymd("20200226"), units="days"))
GathRestrictAnyDays0 <- as.numeric(difftime(ymd(GathRestrictRecomAny), ymd("20200226"), units="days"))
StayAtHomeDays0 <- as.numeric(difftime(ymd(StayAtHome), ymd("20200226"), units="days"))

## Days to plot
DayRange <- 12:60

## Compute cumulative counts of events by type
cumulativeSC <- rep(NA, length(DayRange))
for (i in 1:length(DayRange)) {
  cumulativeSC[i] <- sum(na.omit(BarRestrictDays0)<=DayRange[i])
}

cumulativeRC <- rep(NA, length(DayRange))
for (i in 1:length(DayRange)) {
  cumulativeRC[i] <- sum(na.omit(RestaurantRestrictDays0)<=DayRange[i])
}

cumulativeBC <- rep(NA, length(DayRange))
for (i in 1:length(DayRange)) {
  cumulativeBC[i] <- sum(na.omit(AnyBusinessCloseDays0)<=DayRange[i])
}

cumulativeGR <- rep(NA, length(DayRange))
for (i in 1:length(DayRange)) {
  cumulativeGR[i] <- sum(na.omit(GathRestrictAnyDays0)<=DayRange[i])
}

cumulativeSP <- rep(NA, length(DayRange))
for (i in 1:length(DayRange)) {
  cumulativeSP[i] <- sum(na.omit(StayAtHomeDays0)<=DayRange[i])
}

initCasesSC <- cumulativeSC[length(cumulativeSC)]
initCasesRC <- cumulativeRC[length(cumulativeRC)]
initCasesBC <- cumulativeBC[length(cumulativeBC)]
initCasesGR <- cumulativeGR[length(cumulativeGR)]
initCasesSP <- cumulativeSP[length(cumulativeSP)]

BarRestrictDays0 <- weeksPast(BarRestrictDays0)
RestaurantRestrictDays0 <- weeksPast(RestaurantRestrictDays0)
AnyBusinessCloseDays0 <- weeksPast(AnyBusinessCloseDays0)
GathRestrictAnyDays0 <- weeksPast(GathRestrictAnyDays0)
StayAtHomeDays0 <- weeksPast(StayAtHomeDays0)



## Now recompute to cumulate first easings
## First compute cumulative cases at risk by start of easing phase (April 6 2020)

## Attach data for convenience
attach(StatePolicyStopTableWide)

## Convert to days since 26 February 2020 (first reported US community transmission of COVID-19) 

# here I renamed stayathome and restaurantrestrict (they were previously shelterinplace and restaurantclose)
BarRestrictDays <- as.numeric(difftime(ymd(BarRestrict), ymd("20200226"), units="days"))
RestaurantRestrictDays <- as.numeric(difftime(ymd(RestaurantRestrict), ymd("20200226"), units="days"))
AnyBusinessCloseDays <- as.numeric(difftime(ymd(AnyBusinessClose), ymd("20200226"), units="days"))
GathRestrictAnyDays <- as.numeric(difftime(ymd(GathRestrictRecomAny), ymd("20200226"), units="days"))
StayAtHomeDays <- as.numeric(difftime(ymd(StayAtHome), ymd("20200226"), units="days"))

## Convert to weeks for plotting
BarRestrictDays <- weeksPast(BarRestrictDays)
RestaurantRestrictDays <- weeksPast(RestaurantRestrictDays)
AnyBusinessCloseDays <- weeksPast(AnyBusinessCloseDays)
GathRestrictAnyDays <- weeksPast(GathRestrictAnyDays)
StayAtHomeDays <- weeksPast(StayAtHomeDays)



## Weeks to plot
DayRange <- 4:18


## Compute cumulative counts of events by type
cumulativeSC <- rep(NA, length(DayRange))
for (i in 1:length(DayRange)) {
  cumulativeSC[i] <- sum(na.omit(BarRestrictDays)<=DayRange[i])
}

cumulativeRC <- rep(NA, length(DayRange))
for (i in 1:length(DayRange)) {
  cumulativeRC[i] <- sum(na.omit(RestaurantRestrictDays)<=DayRange[i])
}

cumulativeBC <- rep(NA, length(DayRange))
for (i in 1:length(DayRange)) {
  cumulativeBC[i] <- sum(na.omit(AnyBusinessCloseDays)<=DayRange[i])
}

cumulativeGR <- rep(NA, length(DayRange))
for (i in 1:length(DayRange)) {
  cumulativeGR[i] <- sum(na.omit(GathRestrictAnyDays)<=DayRange[i])
}

cumulativeSP <- rep(NA, length(DayRange))
for (i in 1:length(DayRange)) {
  cumulativeSP[i] <- sum(na.omit(StayAtHomeDays)<=DayRange[i])
}


tileSteps <- function(x,y) {
  x0 <- x[1]
  y0 <- y[1]
  n <- length(x)
  for (i in 2:n) {
    if (y[i]!=y[i-1]) {
      x0 <- c(x0, x[i], x[i])
      y0 <- c(y0, y[i-1], y[i])
    } else {
      x0 <- c(x0, x[i])
      y0 <- c(y0, y[i])
    }
  }
  list(x=x0, y=y0)
}


## Gathering restrictions:  convert cumulative counts to tile traces

##Remove repeated circles
DayRange0 <- DayRange
##DayRange0[c(FALSE, 0==diff(initCasesGR - cumulativeGR))] <- NA

GRgrowth <- lineplot(x=DayRange,
                     ##tileSteps(x=DayRange, y=initCasesGR - cumulativeGR)$x,
                     y=initCasesGR - cumulativeGR,
                     ##tileSteps(x=DayRange, y=initCasesGR - cumulativeGR)$y,
                     lex=1.75,
                     col=niceblue,
                     alpha=1,
                     plot=1
)

GRgrowthPts <- pointsTile(x=DayRange0,
                          y=initCasesGR - cumulativeGR,
                          pch=16,
                          size=circleSize, #1.3
                          alpha=1,
                          col=niceblue,
                          clip="off",
                          plot=1
)

GRgrowthLab <- textTile(labels=initCasesGR - cumulativeGR,
                        x=DayRange0,
                        y=initCasesGR - cumulativeGR,
                        pch=16,
                        col="white",
                        layer=1,
                        alpha=1,
                        cex=textSize,  #0.625
                        plot=1
)

## GR: Add state labels for each day a state adds a restriction
##  (These need to be manually cleaned up in Adobe Illustrator!)
activeDays <- sort(unique(na.omit(GathRestrictAnyDays)))
startY <- rep(20, length(activeDays))

collectedStates <- NULL
StateLabCollectedGR <- StateLineCollectedGR <- vector("list", length(activeDays))
for (i in 1:length(activeDays)) {
  currStates <- na.omit(StatePostal[GathRestrictAnyDays==activeDays[i]])
  collectedStates <- c(collectedStates, currStates)
  currRC <- initCasesGR - cumulativeGR[DayRange==activeDays[i]]
  nStates <- length(currStates)
  currY <- seq(startY[i], by=yOffset, length.out=nStates)

  StateLabCollectedGR[[i]] <- textTile(labels=rev(currStates),
                                       x=rep(activeDays[i]+xOffset, nStates),
                                       y=initCasesGR- currY,
                                       cex=cex,
                                       alpha=1,
                                       col=niceblue,
                                       clip="off",
                                       plot=1)

  if (currY[1]>currRC) {
    startStateLine <- currRC + 2.0
    endStateLine <- currY[1] + 1
  } else {
    startStateLine <- currRC - 2.0
    endStateLine <- currY[length(currY)] - 1
  }
  
  StateLineCollectedGR[[i]] <- linesTile(x=rep(activeDays[i], 2),
                                         y=c(startStateLine, endStateLine),
                                         lex=1.325,
                                         alpha=1,
                                         col=brewer.pal(7,"Blues")[2],
                                         clip="off",
                                         layer=20,
                                         plot=1)
}

## Show uneased, never states
uneasedStates <- FindNotStopStateFunc("GathRestrictRecomAny")
if (length(uneasedStates)==0) uneasedStates <- NA
neverStates <- FindNeverEnactStateFunc("GathRestrictRecomAny")
if (length(neverStates)==0) neverStates <- NA

currY <- seq(45, by=yOffset, length.out=length(uneasedStates))
uneasedStatesGR <- textTile(labels=uneasedStates,
                            x=rep(3.5 + xOffset, length(uneasedStates)),
                            y=currY,
                            cex=cex,
                            alpha=1,
                            col=niceorange,
                            clip="off",
                            plot=1)

currY <- seq(45, by=yOffset, length.out=length(neverStates))
neverStatesGR <- textTile(labels=neverStates,
                          x=rep(5 + xOffset, length(neverStates)),
                          y=currY,
                          cex=cex,
                          alpha=1,
                          col=niceorange,
                          clip="off",
                          plot=1)




## Bar Restriction:  convert cumulative counts to tile traces

##Remove repeated circles
DayRange0 <- DayRange
##DayRange0[c(FALSE, 0==diff(initCasesSC - cumulativeSC))] <- NA

SCgrowth <- lineplot(x=DayRange,
                     ##tileSteps(x=DayRange, y=initCasesSC - cumulativeSC)$x,
                     y=initCasesSC - cumulativeSC,
                     ##tileSteps(x=DayRange, y=initCasesSC - cumulativeSC)$y,
                     lex=1.75,
                     col=niceblue,
                     alpha=1,
                     plot=2
)

SCgrowthPts <- pointsTile(x=DayRange0,
                          y=initCasesSC - cumulativeSC,
                          pch=16,
                          size=circleSize, #1.3
                          alpha=1,
                          col=niceblue,
                          clip="off",
                          plot=2
)

SCgrowthLab <- textTile(labels=initCasesSC - cumulativeSC,
                        x=DayRange0,
                        y=initCasesSC - cumulativeSC,
                        pch=16,
                        col="white",
                        layer=1,
                        alpha=1,
                        cex=textSize,  #0.625
                        plot=2
)

## SC: Add state labels for each day a state adds a restriction
##  (These need to be manually cleaned up in Adobe Illustrator!)
activeDays <- sort(unique(na.omit(BarRestrictDays)))
startY <- rep(20, length(activeDays))

collectedStates <- NULL
StateLabCollectedSC <- StateLineCollectedSC <- vector("list", length(activeDays))
for (i in 1:length(activeDays)) {
  currStates <- na.omit(StatePostal[BarRestrictDays==activeDays[i]])
  collectedStates <- c(collectedStates, currStates)
  currRC <- initCasesSC - cumulativeSC[DayRange==activeDays[i]]
  nStates <- length(currStates)
  currY <- seq(startY[i], by=yOffset, length.out=nStates)
  
  StateLabCollectedSC[[i]] <- textTile(labels=rev(currStates),
                                       x=rep(activeDays[i]+xOffset, nStates),
                                       y=initCasesSC- currY,
                                       cex=cex,
                                       alpha=1,
                                       col=niceblue,
                                       clip="off",
                                       plot=2)
  
  if (currY[1]>currRC) {
    startStateLine <- currRC + 2.0
    endStateLine <- currY[1] + 1
  } else {
    startStateLine <- currRC - 2.0
    endStateLine <- currY[length(currY)] - 1
  }
  
  StateLineCollectedSC[[i]] <- linesTile(x=rep(activeDays[i], 2),
                                         y=c(startStateLine, endStateLine),
                                         lex=1.325,
                                         alpha=1,
                                         col=brewer.pal(7,"Blues")[2],
                                         clip="off",
                                         layer=20,
                                         plot=2)
}

## Show uneased, never states
uneasedStates <- FindNotStopStateFunc("BarRestrict")
if (length(uneasedStates)==0) uneasedStates <- NA
neverStates <- FindNeverEnactStateFunc("BarRestrict")
if (length(neverStates)==0) neverStates <- NA

currY <- seq(45, by=yOffset, length.out=length(uneasedStates))
uneasedStatesSC <- textTile(labels=uneasedStates,
                            x=rep(3.5 + xOffset, length(uneasedStates)),
                            y=currY,
                            cex=cex,
                            alpha=1,
                            col=niceorange,
                            clip="off",
                            plot=2)

currY <- seq(45, by=yOffset, length.out=length(neverStates))
neverStatesSC <- textTile(labels=neverStates,
                          x=rep(5 + xOffset, length(neverStates)),
                          y=currY,
                          cex=cex,
                          alpha=1,
                          col=niceorange,
                          clip="off",
                          plot=2)




## Restaurant restrictions:  convert cumulative counts to tile traces

##Remove repeated circles
DayRange0 <- DayRange
##DayRange0[c(FALSE, 0==diff(initCasesRC - cumulativeRC))] <- NA

RCgrowth <- lineplot(x=DayRange,
                     ##tileSteps(x=DayRange, y=initCasesRC - cumulativeRC)$x,
                     y=initCasesRC - cumulativeRC,
                     ##tileSteps(x=DayRange, y=initCasesRC - cumulativeRC)$y,
                     lex=1.75,
                     col=niceblue,
                     alpha=1,
                     plot=3
)

RCgrowthPts <- pointsTile(x=DayRange0,
                          y=initCasesRC - cumulativeRC,
                          pch=16,
                          size=circleSize, #1.3
                          alpha=1,
                          col=niceblue,
                          clip="off",
                          plot=3
)

RCgrowthLab <- textTile(labels=initCasesRC - cumulativeRC,
                        x=DayRange0,
                        y=initCasesRC - cumulativeRC,
                        pch=16,
                        col="white",
                        layer=1,
                        alpha=1,
                        cex=textSize,  #0.625
                        plot=3
)

## RC: Add state labels for each day a state adds a restriction
##  (These need to be manually cleaned up in Adobe Illustrator!)
activeDays <- sort(unique(na.omit(RestaurantRestrictDays)))
startY <- rep(20, length(activeDays))

collectedStates <- NULL
StateLabCollectedRC <- StateLineCollectedRC <- vector("list", length(activeDays))
for (i in 1:length(activeDays)) {
  currStates <- na.omit(StatePostal[RestaurantRestrictDays==activeDays[i]])
  collectedStates <- c(collectedStates, currStates)
  currRC <- initCasesRC - cumulativeRC[DayRange==activeDays[i]]
  nStates <- length(currStates)
  currY <- seq(startY[i], by=yOffset, length.out=nStates)
  
  StateLabCollectedRC[[i]] <- textTile(labels=rev(currStates),
                                       x=rep(activeDays[i]+xOffset, nStates),
                                       y=initCasesRC- currY,
                                       cex=cex,
                                       alpha=1,
                                       col=niceblue,
                                       clip="off",
                                       plot=3)
  
  if (currY[1]>currRC) {
    startStateLine <- currRC + 2.0
    endStateLine <- currY[1] + 1
  } else {
    startStateLine <- currRC - 2.0
    endStateLine <- currY[length(currY)] - 1
  }
  
  StateLineCollectedRC[[i]] <- linesTile(x=rep(activeDays[i], 2),
                                         y=c(startStateLine, endStateLine),
                                         lex=1.325,
                                         alpha=1,
                                         col=brewer.pal(7,"Blues")[2],
                                         clip="off",
                                         layer=20,
                                         plot=3)
}

## Show uneased, never states
uneasedStates <- FindNotStopStateFunc("RestaurantRestrict")
if (length(uneasedStates)==0) uneasedStates <- NA
neverStates <- FindNeverEnactStateFunc("RestaurantRestrict")
if (length(neverStates)==0) neverStates <- NA

currY <- seq(45, by=yOffset, length.out=length(uneasedStates))
uneasedStatesRC <- textTile(labels=uneasedStates,
                            x=rep(5 + xOffset, length(uneasedStates)),
                            y=currY,
                            cex=cex,
                            alpha=1,
                            col=niceorange,
                            clip="off",
                            plot=3)

currY <- seq(45, by=yOffset, length.out=length(neverStates))
neverStatesRC <- textTile(labels=neverStates,
                          x=rep(7 + xOffset, length(neverStates)),
                          y=currY,
                          cex=cex,
                          alpha=1,
                          col=niceorange,
                          clip="off",
                          plot=3)




## Nonessential business closures:  convert cumulative counts to tile traces

##Remove repeated circles
DayRange0 <- DayRange
##DayRange0[c(FALSE, 0==diff(initCasesBC - cumulativeBC))] <- NA

BCgrowth <- lineplot(x=DayRange,
                     ##tileSteps(x=DayRange, y=initCasesBC - cumulativeBC)$x,
                     y=initCasesBC - cumulativeBC,
                     ##tileSteps(x=DayRange, y=initCasesBC - cumulativeBC)$y,
                     lex=1.75,
                     col=niceblue,
                     alpha=1,
                     plot=4
)

BCgrowthPts <- pointsTile(x=DayRange0,
                          y=initCasesBC - cumulativeBC,
                          pch=16,
                          size=circleSize, #1.3
                          alpha=1,
                          col=niceblue,
                          clip="off",
                          plot=4
)

BCgrowthLab <- textTile(labels=initCasesBC - cumulativeBC,
                        x=DayRange0,
                        y=initCasesBC - cumulativeBC,
                        pch=16,
                        col="white",
                        layer=1,
                        alpha=1,
                        cex=textSize,  #0.625
                        plot=4
)

## BC: Add state labels for each day a state adds a restriction
##  (These need to be manually cleaned up in Adobe Illustrator!)
activeDays <- sort(unique(na.omit(AnyBusinessCloseDays)))
startY <- rep(20, length(activeDays))

collectedStates <- NULL
StateLabCollectedBC <- StateLineCollectedBC <- vector("list", length(activeDays))
for (i in 1:length(activeDays)) {
  currStates <- na.omit(StatePostal[AnyBusinessCloseDays==activeDays[i]])
  collectedStates <- c(collectedStates, currStates)
  currRC <- initCasesBC - cumulativeBC[DayRange==activeDays[i]]
  nStates <- length(currStates)
  currY <- seq(startY[i], by=yOffset, length.out=nStates)
  
  StateLabCollectedBC[[i]] <- textTile(labels=rev(currStates),
                                       x=rep(activeDays[i]+xOffset, nStates),
                                       y=initCasesBC- currY,
                                       cex=cex,
                                       alpha=1,
                                       col=niceblue,
                                       clip="off",
                                       plot=4)
  
  if (currY[1]>currRC) {
    startStateLine <- currRC + 2.0
    endStateLine <- currY[1] + 1
  } else {
    startStateLine <- currRC - 2.0
    endStateLine <- currY[length(currY)] - 1
  }
  
  StateLineCollectedBC[[i]] <- linesTile(x=rep(activeDays[i], 2),
                                         y=c(startStateLine, endStateLine),
                                         lex=1.325,
                                         alpha=1,
                                         col=brewer.pal(7,"Blues")[2],
                                         clip="off",
                                         layer=20,
                                         plot=4)
}

## Show uneased, never states
uneasedStates <- FindNotStopStateFunc("AnyBusinessClose")
if (length(uneasedStates)==0) uneasedStates <- NA
neverStates <- FindNeverEnactStateFunc("AnyBusinessClose")
if (length(neverStates)==0) neverStates <- NA

currY <- seq(45, by=yOffset, length.out=length(uneasedStates))
uneasedStatesBC <- textTile(labels=uneasedStates,
                            x=rep(3.5 + xOffset, length(uneasedStates)),
                            y=currY,
                            cex=cex,
                            alpha=1,
                            col=niceorange,
                            clip="off",
                            plot=4)

currY <- seq(45, by=yOffset, length.out=length(neverStates))
neverStatesBC <- textTile(labels=neverStates,
                          x=rep(5 + xOffset, length(neverStates)),
                          y=currY,
                          cex=cex,
                          alpha=1,
                          col=niceorange,
                          clip="off",
                          plot=4)



## Stay at home:  convert cumulative counts to tile traces

##Remove repeated circles
DayRange0 <- DayRange
##DayRange0[c(FALSE, 0==diff(initCasesSP - cumulativeSP))] <- NA

SPgrowth <- lineplot(x=DayRange,
                     ##tileSteps(x=DayRange, y=initCasesSP - cumulativeSP)$x,
                     y=initCasesSP - cumulativeSP,
                     ##tileSteps(x=DayRange, y=initCasesSP - cumulativeSP)$y,
                     lex=1.75,
                     col=niceblue,
                     alpha=1,
                     plot=5
)

SPgrowthPts <- pointsTile(x=DayRange0,
                          y=initCasesSP - cumulativeSP,
                          pch=16,
                          size=circleSize, #1.3
                          alpha=1,
                          col=niceblue,
                          clip="off",
                          plot=5
)

SPgrowthLab <- textTile(labels=initCasesSP - cumulativeSP,
                        x=DayRange0,
                        y=initCasesSP - cumulativeSP,
                        pch=16,
                        col="white",
                        layer=1,
                        alpha=1,
                        cex=textSize,  #0.625
                        plot=5
)

## SP: Add state labels for each day a state adds a restriction
##  (These need to be manually cleaned up in Adobe Illustrator!)
activeDays <- sort(unique(na.omit(StayAtHomeDays)))
startY <- rep(20, length(activeDays))

collectedStates <- NULL
StateLabCollectedSP <- StateLineCollectedSP <- vector("list", length(activeDays))
for (i in 1:length(activeDays)) {
  currStates <- na.omit(StatePostal[StayAtHomeDays==activeDays[i]])
  collectedStates <- c(collectedStates, currStates)
  currRC <- initCasesSP - cumulativeSP[DayRange==activeDays[i]]
  nStates <- length(currStates)
  currY <- seq(startY[i], by=yOffset, length.out=nStates)
  
  StateLabCollectedSP[[i]] <- textTile(labels=rev(currStates),
                                       x=rep(activeDays[i]+xOffset, nStates),
                                       y=initCasesSP- currY,
                                       cex=cex,
                                       alpha=1,
                                       col=niceblue,
                                       clip="off",
                                       plot=5)
  
  if (currY[1]>currRC) {
    startStateLine <- currRC + 2.0
    endStateLine <- currY[1] + 1
  } else {
    startStateLine <- currRC - 2.0
    endStateLine <- currY[length(currY)] - 1
  }
  
  StateLineCollectedSP[[i]] <- linesTile(x=rep(activeDays[i], 2),
                                         y=c(startStateLine, endStateLine),
                                         lex=1.325,
                                         alpha=1,
                                         col=brewer.pal(7,"Blues")[2],
                                         clip="off",
                                         layer=20,
                                         plot=5)
}

## Show uneased, never states
uneasedStates <- FindNotStopStateFunc("StayAtHome")
if (length(uneasedStates)==0) uneasedStates <- NA
neverStates <- FindNeverEnactStateFunc("StayAtHome")
if (length(neverStates)==0) neverStates <- NA

currY <- seq(45, by=yOffset, length.out=length(uneasedStates))
uneasedStatesSP <- textTile(labels=uneasedStates,
                            x=rep(3.5 + xOffset, length(uneasedStates)),
                            y=currY,
                            cex=cex,
                            alpha=1,
                            col=niceorange,
                            clip="off",
                            plot=5)

currY <- seq(45, by=yOffset, length.out=length(neverStates))
neverStatesSP <- textTile(labels=neverStates,
                          x=rep(5 + xOffset, length(neverStates)),
                          y=currY,
                          cex=cex,
                          alpha=1,
                          col=niceorange,
                          clip="off",
                          plot=5)




##### cumulative scores plot

## Check whether ever implemented
allMeas0 <- cbind(BarRestrictDays0, RestaurantRestrictDays0,
                  AnyBusinessCloseDays0, GathRestrictAnyDays0,
                  StayAtHomeDays0)

## Check when first eased 
allMeas <- cbind(BarRestrictDays, RestaurantRestrictDays,
                 AnyBusinessCloseDays, GathRestrictAnyDays,
                 StayAtHomeDays)

## Replace never implemented cases with -10000
for (i in 1:nrow(allMeas)) {
  for (j in 1:ncol(allMeas)) {
    if (is.na(allMeas0[i,j]))
      allMeas[i,j] <- -10000
  }
}

## Replace never eased cases with 10000
allMeas[is.na(allMeas)] <- 10000

res <- matrix(NA, nrow=length(DayRange), ncol=6)
for (i in 1:length(DayRange)) {
  stateCounts <- apply(allMeas, 1, function(x){sum(x<=DayRange[i])})
  ##stateCounts <- abs(stateCounts-5)  ## convert to "measures lacking" to make plotting easier
  for (j in 0:5) {
    res[i,(j+1)] <- sum(stateCounts==j)
  }
}


ctr <- 0
polycount <- sum(res>0)

ycs <- t(apply(cbind(rep(0,length(DayRange)),res),1,cumsum))
NmeasCollectedTrace <- vector("list", polycount)
reds <- rev(brewer.pal(6,"Blues")[1:6])
for (i in 1:length(DayRange)) {
  for (j in 1:6) {
    xvert <- c(DayRange[i]-0.5, DayRange[i]+0.5,
               DayRange[i]+0.5, DayRange[i]-0.5)
    yvert <- c(ycs[i,j], ycs[i,j], ycs[i,j+1], ycs[i,j+1])
    if (ycs[i,j+1]> ycs[i,j]) {
      ctr <- ctr + 1
      NmeasCollectedTrace[[ctr]] <- polygonTile(x=xvert,
                                                y=yvert,
                                                col=NA,
                                                fill=reds[j],
                                                layer=8,
                                                alpha=1,
                                                plot=6
      )
    }
  }
}

cumulativeLabelTrace <- textTile(labels=c("All eased\n or absent", "One\n uneased", "Two\n uneased", "Three\n uneased", "Four\n uneased", "All five\n uneased"),
                                 x=c(15, 15, 14.5, 6, 6, 6),
                                 y=c(44, 30, 10, 45, 30, 10),
                                 cex=c(rep(1.0,5), 1.0),
                                 col=c(rep(reds[1], 3), rep(reds[6], 3)),
                                 plot=6
)

##

datestamp <- textTile(labels=paste("\uA9 COVID-19 State Policy Team / Christopher Adolph & Erika Steiskal     ",
                                   format(Sys.time(), "%a %d %b %Y %H:%M:%S"),
                                   "    http://covid19statepolicy.org"),
                      x=10,
                      y=-25,
                      cex=0.5,
                      clip="off",
                      plot=5)


## Send everything to tile for plotting
##file <- "Figure1prepolish"
if (fancy) {
  output <- list(file=file, width=10, family="GillSans") 
} else {
  output <- list(file=file, width=10) 
}
tile(GRgrowth, GRgrowthPts, GRgrowthLab,
     StateLabCollectedGR, StateLineCollectedGR,
     SCgrowth, SCgrowthPts, SCgrowthLab,
     StateLabCollectedSC, StateLineCollectedSC,
     RCgrowth, RCgrowthPts, RCgrowthLab,
     StateLabCollectedRC, StateLineCollectedRC,
     BCgrowth, BCgrowthPts, BCgrowthLab,
     StateLabCollectedBC, StateLineCollectedBC,
     SPgrowth, SPgrowthPts, SPgrowthLab,
     StateLabCollectedSP, StateLineCollectedSP,
     neverStatesGR, neverStatesSC, neverStatesRC,
     neverStatesBC, neverStatesSP,
     uneasedStatesGR, uneasedStatesSC, uneasedStatesRC,
     uneasedStatesBC, uneasedStatesSP,
     NmeasCollectedTrace, cumulativeLabelTrace,
     datestamp,
     RxC=c(3,2),
     limits=c(min(DayRange) - 1, max(DayRange) + 1, -3, 51), 
     xaxis=list(at=seq(2,18,2), y=0,
                labels=c("Weeks since 1 April:", "", "6", "8", "10", "12","14", "16", "18")), #seq(60,120,30)
     yaxis=list(log=FALSE, major=FALSE,
                at=c(0, 10, 20, 30, 40, 50)),
     yaxistitle=list(labels1="States Maintaining Uneased Policy",
                     labels3="States Maintaining Uneased Policy",
                     labels5="States Maintaining Uneased Policy", x=0.0),
     plottitle=list(labels1="GATHERING RESTRICTIONS",
                    labels2="BAR RESTRICTIONS",
                    labels3="RESTAURANT RESTRICTIONS",
                    labels4="ANY BUSINESS CLOSURES",
                    labels5="STAY-AT-HOME ORDERS",
                    labels6="CUMULATIVE MEASURES UNEASED",
                    y=0.9),
     xaxistitle=list(labels=" "), ##Weeks Since 1 March 2020
     height=list(xaxistitle=3, plottile=3, spacer=2.75, plot="golden"), 
     width=list(yaxis.labelspace=-0.15),
     output=output
)
if (fancy) embedGS(paste(file,".pdf",sep=""))

## Then import to Illustrator for many manual polishing steps to make readable
## Including manual adjustment of the states in orange based on details from the underlying COVID-19 state policy data

